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ABSTRACT 

We report the radial-velocity discovery of a second planetary mass companion to the K0 V star 
HD 37605, which was already known to host an eccentric, P ~ 55 days Jovian planet, HD 376056. 
This second planet, HD 37605c, has a period of ~ 7.5 years with a low eccentricity and an Msini of 
~ 3.4 Afj up . Our discovery was made with the nearly 8 years of radial velocity follow-up at the Hobby- 
Eberly Telescope and Keck Observatory, including observations made as part of the Transit Ephemeris 
Refinement and Monitoring Survey (TERMS) effort to provide precise ephemerides to long-period 
planets for transit follow-up. With a total of 137 radial velocity observations covering almost eight 
years, we provide a good orbital solution of the HD 37605 system, and a precise transit ephemeris for 
HD 376056. Our dynamic analysis reveals very minimal planet-planet interaction and an insignificant 
transit time variation. Using the predicted ephemeris, we performed a transit search for HD 376056 
with the photometric data taken by the T12 0.8-m Automatic Photoelectric Telescope (APT) and the 
Microvariability and Oscillations of Stars (MOST) satellite. Though the APT photometry did not 
capture the transit window, it characterized the stellar activity of HD 37605, which is consistent of 
it being an old, inactive star, with a tentative rotation period of 57.67 days. The MOST photometry 
enabled us to report a dispositive null detection of a non- grazing transit for this planet. Within 
the predicted transit window, we exclude an edge-on predicted depth of 1.9% at the 3> 10cr level, 
and exclude any transit with an impact parameter 6 > 0.951 at greater than 5er. We present the 
BOOTTRAN package for calculating Keplerian orbital parameter uncertainties via bootstrapping. 
We made a comparison and found consistency between our orbital fit parameters calculated by the 
RVLIN package and error bars by BOOTTRAN with those produced by a Bayesian analysis using 
MCMC. 

Subject headings: planetary systems — stars: individual (HD 37605) 
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1. INTRODUCTION 
1.1. Context 

Jupiter analogs orbiting other stars represent the first 
signposts of true Solar System analogs, and the eccentric- 
ity distribution of these planets with a > 3 AU will reveal 
how rare or frequent true Jupiter analogs are. To date, 
only 9 "Jupiter analogs" have been well-characterized in 
the peer reviewed literature[3 (defined he re as P > 8 
years , 4 > Msini > 0.5 M Jup , and e < 0.3; lWright et al.l 
120111 exoplanets.org). As the duration of existing planet 
searches approach 10-20 years, more and more Jupiter 
analogs will emerge from their longest-obs erved targets 
([Wittenrnver et al.|[20ll IBoisse et al.ll2012t) . 

Of the over 700 exoplanets discovered to date, 
nearly 200 are kn own to transit t heir host star 
(Wri ght et al.l 120111 exoplanets.org; iSchneider et al.l 
120111 exoplanet.eul and many thousands more candi- 
dates have been discovered by the Kepler telescope. Of 
all of these planets, only three orbit stars with V < 8F^I 
and all have P < 4 days. Long period planets are less 
likely than close-in planets to transit unless their orbits 
are highly eccentric and favorably oriented, and indeed 
only 2 transiting planets with P > 20 days have been 
discovered around st ars with V < 10, and both have 
e > .65 (HD 8 0606. ILaughlin et al I [20091 iFossev et all 
[2009t HD 17156. lFrscher et alll2007l iBarbieri et alll2007l: 
both highly eccentric systems were discovered first with 
radial velocities) . 

Long period planets not known to transit can have 
long transit windows due to both the large duration of 
any edge-on transit and higher phase uncertainties (since 
such uncertainties scale with the period of the orbit). 
Long term radial velocity monitoring of stars, for in- 
stance for the discovery of low amplitude signals, can 
produce collateral benefits in the form of orbit refine- 
ment for a tr ansit search and the identification of Jupiter 
analogs fe.g,. IWright et al1l2009D . Herein, we describe an 
example of both. 

1.2. Initial Discovery and Followup 

The inner planet in the system, HD 376056, was the 
first planet discovered with the Ho bby-Eberly Telescop e 
(HET) at McDonald Observatory (|Cochran et al.ll2004D . 
It is a super Jupiter (Msini = 2.41 Mj up ) on an eccen- 
tric orbit e = 0.67 wi th an orbital period in the "period 
valley" (P = 55 davs; IWright et alj|2009h . 

W.C., M.E., and P.J.M. of the University of Texas at 
Austin, continued observations in order to get a much 
better orbit determination and to begin searching for 
transits. With the first new data in the fall of 2004, 
it became obvious that another perturber was present in 
the system, first from a trend in the r adial velocity (RV) 
resid uals (i.e., a none-zero dv/dt; IWittenmver et all 
120071 ). and later from curvature in the residuals. By 
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2009, the residuals to a one-planet fit were giving rea- 
sonable constraints on the orbit of a second planet, HD 
37605c, and by early 2011 the orbital parameters of the c 
component were clear, and the Texas team was preparing 
the system for publication. 

1.3. TERMS Data 

The Transit Ephem eris Refinement and Monitor- 
ing Survey (TERMS; [Kane et al.l I2009D seeks to re- 
fine the ephemerides of the known exoplanets orbit- 
ing bright, nearby stars with sufficient precision to ef- 
ficiently search for the planetary transits of planets 
with periast ron distances greater than a few hundredths 
of an AU (Kane et all 1201 Hi iPilvavskv et"all I2011ah 
iDragomir et al.ll2011l ). This will provide the radii of plan- 
ets not experiencing continuous high levels of insolation 
around nearby, easily studied stars. 

In 2010, S.M. and J.T.W. began radial velocity obser- 
vations of HD 376056 at HET from Penn State Univer- 
sity for TERMS, to refine the orbit of that planet for 
a future transit search. These observations, combined 
with Keck radial velocities from the California Planet 
Survey (CPS) consortium from 2006 onward, revealed 
that there was substantia l curvature to the radi al velocity 
residuals to the original ICochran et all (|2004D solution. 
In October 2010 monitoring was intensified at HET and 
at Keck Observatory by A.W.H., G.W.M., J.T.W., and 
H.I., and with these new RV data and the pre v iously 
published measurements from IWittenmver et all (|2007| ) 
they obtained a preliminary solution for the outer planet. 
The discrepancy between the original orbital fit and the 
new fit (assuming one planet) was presented at the Jan- 
uary 2011 meeting of the American Astronomical Society 
(jKane et all(2^Tl3) . 

1.4. Synthesis and Outline 

In early 2011, the Texas and TERMS teams combined 
efforts and began joint radial velocity analysis, dynam- 
ical modeling , spectroscopic a nalysis, and photometric 
observations (jKane et al.ll2012j ). The resulting complete 
two-planet orbital solution allows for a sufficiently pre- 
cise transit ephemeris for the 6 component to be calcu- 
lated for a thorough transit search. We herein report the 
transit exclusion of HD 376056 and a stable dynamical 
solution to the system. 

In § [21 we describe our spectroscopic observations and 
analysis, which provided the radial velocities and the stel- 
lar properties of HD 37605. § [3] details the orbital solu- 
tion for the HD 37605 system, including a comparison 
with MCMC Keplerian fits, and our dynamical analysis. 
We report our photometric observations on HD 37605 
and the dispositive null detection^ of non-grazing tran- 
sits of HD 376056 in § H After § EJ Summary and Con- 
clusion, we present updates on M sin i of two previously 
published systems (HD 114762 and HD 168443) in § [EJ 
In the Appendix we describe the algorithm used in the 
package BOOTTRAN (for calculating orbital parameter 
error bars; see § 13.21). 



14 55 CncejMcArtlmr ct al. 2004,; IDemorv et "all |20TlD , HD 
189733 (TBouchv et alJ 120051). and HD 209458 I IHenrv et al.1 feoOOl ; 



Charbonncau ct al. 2000). 



15 A dispositive null detection is one that disposes of the ques- 
tion of whether an effect is present, as opposed to one that merely 
fails to detect a purported or hypothetical effect that may yet lie 
beneath the detection threshold. The paragon of dispositive null 
detections is the Michelson -Morley demonstration that the luminif- 
erous ether does not exist (Michclson & Morlcy 1887). 
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2. SPECTROSCOPIC OBSERVATIONS AND ANALYSIS 
2.1. HET and Keck Observations 

Observations on HD 37605 at HET started Decem- 
ber of 2003. In total, 101 RV observations took place 
over the course of almost eight years, taking advantage 
of the queue scheduling capabilities of HET. The queue 
scheduling of HET allows for small amounts of telescope 
time to be optimally used throughout the year, and for 
new observing priorities to be implemented immediately, 
rather than on nex t allocated night or after TAC and 
scheduling process (jShetrone et al J 12007). The observa- 
tions were t aken throu gh the High Resolution Spectro- 
graph (HRS; lTulllll998l ) situated at the basement of the 
HET building. This fiber-fed spectrograp h has a typica l 
long-term Doppler error of 3 - 5 m s" 1 (IBaluevI 120091) . 
The observations were taken with the spectrograph con- 
figured at a resolving power of R =60,000. For more 
details, see iCochran et all (pOOl . 

Observations at Keck were taken starting August 
2006. A set of 33 observations spanning ove r five years 
were made through the HIRES spectrometer (|Vogt et al.l 
1994) on the Keck I telescope, which has a long-term 
Dopp ler error of 0.9 - 1.5 m s _1 (e.g. iHoward et al.l 
2009). The observations were taken at a resolving power 
of R =55,000. For more details, see IHoward et~aLl (|2009h 
and lValenti et al.1 (|2009h . 

Both our HET and Keck spectroscopic observations 
were taken with an iodine cell placed in the light path 
to provide wavelength standard and information on the 
instrument response function^ (IRF ) for radial velo city 
extraction (|Marcv fc Butlerlll992tlButler et al.lll996l ). In 
addition, we also have observations taken without iodine 
cell to produce stellar spectrum templates - on HET and 
Keck, respectively. The stellar spectrum templates, after 
being deconvolved with the IRF, are necessary for both 
radial velocity extraction and stellar property analysis. 
The typical working wavelength range for this technique 
is roughly 5000 A- 6000 A. 

2.2. Data Reduction and Doppler Analysis 

In this section, we describe our data reduction and 
Doppler analysis of the HET observations. We re- 
duced the Keck data with th e standard CPS pip eline, 
as described in, for e xample, IHoward et al.l (|2011l) and 
Uohnson et all (|2011bl ). 

We have constructed a complete pipeline for analyz- 
ing HET data - from raw data reduction to radial ve- 
locity extraction. Th e raw reduction is d o ne us ing the 
REDUCE package bv Fiskunov fc Valentil pOOl . This 
package is design ed to optima lly extract echelle spectra 
from 2-D images (|Hornelll986f) . Our pipeline corrects for 
cosmic rays and scattered light. In order to make the 
data reduction process completely automatic, we have 
developed our own algorithm for tracing the echelle or- 
ders of HRS and replaced the original semi-automatic 
algorithm from the REDUCE package. 

After the raw data reduction, the stellar spectrum tem- 
plate is deconvolved using IRF derived from an iodine 
flat on the night of observation. There were two decon- 
volved stellar spectrum templates (DSST) derived from 

16 Some authors refer to this as the "point spread function" or 
the "instrumental profile" of the spectrograph. 



TABLE 1 
STELLAR PARAMETERS 



Parameter 


Value 


Spectral typc a 


K0 V 




44.0 i 2.1 


V 


8.661 ± 0.013 


T eS (K) 


5448 ±44 


logs 


4.511 ± 0.024 


[Fe/H] 


0.336 ± 0.030 


BC 


-0.144 


M bol 


5.301 




0.590 ± 0.058 


R* (Re) 


0.901 ± 0.045 c 


hh (M e ) 


1.000 ± 0.050 c 


v sini 


< 1 km s — 1 


Age b 


~ 7 Gyr 


a IESAI U997D:lvan Leeuwenl j200Sl). 
b llsaacson & Fischer) 120101) , see 5 [4~T1 

c 5% relative errors, not the SME intrinsic errors. 
See footnote fTT] for details. 



HET/HRS observations and one from Keck/HIRES. 
Throughout this work, we use the Keck DSST, which 
is of better quality thanks to a better known IRF of 
HIRES and a superior deconvolution a lgorithm in the 
CPS pipeline (IHoward et al.lf2009l [20TlT> . 

Then the pipeline proceeds with barycentric correc- 
tion and radial velocity extraction for each observa- 
tion. We have adop t ed the Doppler code from CP S (e.g. 
IHoward et all I2009L 12011 Uohnson etaU I2011bl) . The 
code is tailored to be fully functional with HET/HRS- 
formatted spectra, and it is capable of working with ei- 
ther an HET DSST or a Keck one. 

The 101 HET RV observations include 44 observations 
which produ ced the published v e locitie s in lCochran et al.l 
(|2004f) and iWittenmver et al.l (|2007l ). 34 observations 
also done by the Texas team in follow-up work after 2007, 
and 23 observations taken as part of TERMS program. 
We have performed re-reduction on these 44 observations 
together with all the rest 57 HET observations through 
our pipeline. This has the advantage of eliminating one 
free parameter in the Keplerian fit - the offset between 
two Doppler pipelines. 

Two out of the 101 HET observations were excluded 
due to very low average signal-to-noise ratio per pixel 
(< 20), and one observation taken at twilight was also 
rejected as such observation normally results in low accu- 
racy due to the significant contamination by the residual 
solar spectrum (indeed this velocity has a residual of over 
100 m s _1 against best Keplerian fit, much larger than 
the ~ 8 m s _1 RV error). 

All the HET and Keck radial velocities used in this 
work (98 from HET and 33 from Keck) are listed in Table 

El 

2.3. Stellar Analysis 

HD 37605 is a K0 V star (V - 8.7) with high proper 
motion at a dista nce of 44.0 ± 2.1 pc (|ESAI Il997t 
Ivan Leeiiw cn 2008). We derived its stellar properties 
based on analysis on a high-resolution spectrum taken 
with Keck HIRES (without iodine cell in the light path) . 
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Table [T] lists the results of our analysisQ, including the ef- 
fective temperature T c g, surface gravity logg, iron abun- 
dance [Fe/H], projected rotational velocity vsmi, bolo- 
metric correction BC, bolometric magnitude M\, a \, stellar 
luminosity L ± , stellar radius R+, stellar mass and age. 
HD 37605 is found to be a metal rich star ([Fe/H] - 0.34) 
with M* ~ 1.0 M Q and R* ~ O.9i? . 

We followed the procedure d escribed in 
IValenti fc Fischer] ()2005|) and also in iValenti et al.l 
(2009) with improvements. Briefly, the observed 
spectrum is fitted with a synt hetic spectrum using Spec- 
troscopy Made Easy (SME; IValenti fc Piskunovl [l99ll 
to derive T e s, logg, [Fe/H], vs'mi, and so on, which 
are used to derive the bolometric correction BC and 
consequently. Then an isochro ne fit by interpolating 
tabulated Yonsei-Yale isochrones (jDemarque et al.ll2004f) 
using derived stellar parameters from SME is performed 
to calculate M± and loggi so values (along with age and 
stellar radius). Next, Valenti et al. (2009) introduced 
an outside loop which re-runs SME with logg fixed at 
loggi S o, followed by another isochrone fit deriving a new 
log loggi S o using the updated SME results. The loop 
continues until logg values converge. This additional 
iterative procedure to enforce self-consistency on log g is 
shown to im prove the accuracy of other derived stellar 
parameters ()Valenti et al.l 120091 ) . The stellar radius 
and log g reported here in Table [1] are derived from the 
final isochrone fit, which are consistent with the purely 
spectroscopic results. The gravity (logg = 4.51) is also 
consistent with the purely spectroscopic gravity (4.44) 
based on strong Mg b damping wings, so for HD 37605 
t he iteration process is optional. 

ICochran et al.l (|2004[) reported the values of T c ff , logg, 
and [Fc/H] for HD 37605, a nd their estimates ag ree with 
ours within la uncertainty. Sa ntos et al.l (|2005l ) also es- 
timated T c ff, logg, [Fe/H], and M*, all of which agree 
with our values within la. Our stellar mass and radius 
estimates are also consis tent with the ones d erived from 
the empirical method bv lTorres et al.l (|2010l ). 

Our SME analysis indicates that the rotation of the 
star (vsmi) is likely < 1 km s _1 (corresponding to rota- 
tion period > 46 days) . We have used various methods to 
estimate stellar parameters from the spectrum, including 
the incorporation of color and absolute magnitude infor- 
mation and the Mg b triplet to constrain logg, and var- 
ious macroturbulent velocity prescriptions. All of these 
approaches yield results consistent with an undetectable 
level of rotational broadening, with an upper limit of 1-2 
km s , consistent with the tentative photometric period 
57.67 days derived from the APT data (See gH|). 

3. ORBITAL SOLUTION 

3.1. Transit Ephemeris 

The traditional parameters for reporting the 
ephemerides of spectroscopic binaries are P,K,e,uj, 
and T P , the last being the time of periastron passage 
([Wright fc Howard 2009]) . This information is sufficient 

17 Note that the errors on the stellar radius i?* and mass 
listed in Tablc[T]are not intrinsic to the SME code, but are 5%xR+ 
and 5%xM,. This is because the intrinsic errors reported by SME 
do not include the errors stemming from the adopted stellar mod- 
els, and a more realistic precision for and M* would be around 
~ 5%. Intrinsic errors reported by SME are 0.015 Lq for R* and 
0.017 Mq for M*. 



to predict the phase of a planet at any point in the 
future in principle, but the uncertainties in those param- 
eters alone are insufficient to compute the uncertainty 
in orbital phase without detailed knowledge of the 
covariances among the parameters. 

This problem is particularly acute when determining 
transit or secondary eclipse times for planets with near 
circular orbits, where ax p and a^ can be highly covariant. 
In such cases the circular case is often not excluded by 
the data, and so the estimation of e includes the case 
e = 0, where w is undefined. If the best or most likely 
value of e in this case is small but not zero, then it is 
associated with some nominal value of u, but a u will 
be very large (approaching 7r). Since T p represents the 
epoch at which the true anomaly equals 0, T p will have a 
similarly large uncertainty (approaching P), despite the 
fact that the phase of the system may actually be quite 
precisely known! 

In practice even the ephemerides of planets with well 
measured eccentricities suffer from lack of knowledge of 
the covariance in parameters, in particular T p and P 
(whose covariance is sensitive to the approximate epoch 
chosen for T p ). To make matters worse, the nature of 
"ler" uncertainties in the literature is inconsistent. Some 
authors may report uncertainties generated while holding 
all or some other parameters constant (for instance, by 
seeing at what excursion from the nominal value \ 2 is re ~ 
duced by 1), while others using bootstrapping or MCMC 
techniques may report the variance in a parameter over 
the full distribution of trials. In any case, covariances are 
rarely reported, and in some cases authors even report 
the most likely values on a parameter-by-parameter basis 
rather than a representative "best fit" , resulting in a set 
of parameters that is not self-consistent. 

The TERMS strategy for refining ephemerides there- 
fore begins with the recalculation of transit time uncer- 
tainties directly from the archival radial velocity data. 
We used bootstrapping (see Appendix) with the time of 
conjunction, T c (equivalent to transit center, in the case 
of transiting planets) computed independently for each 
trial. For systems whose transit time uncertainty makes 
definitive observations implausible or impossible due to 
the accumulation of errors in phase with time, we sought 
additional RV measurements to "lock down" the phase 
of the planet. 

3.2. The 37605 System 

There are in total 137 radial velocities used in the Kep- 
lerian fit for the HD 37605 system. In addition to the 98 
HET velocities and 33 Keck ones (see N2.2p . we also in- 
cluded siJ^I velocities from ICochran et alj (|2004j ) which 
were derived from observations taken with the McDon- 
ald Observatory 2.1 m Telescope (hereafter the 2.1 m 
telescope). 

We used the RVLIN package by IWright &: Howard! 
(2009) to perform the Keplerian fit. This package is 
based on the Levenberg-Marquardt algorithm and is 
made efficient in searching parameter space by exploiting 
the linear parameters. The uncertainties of the parame- 

18 The velocity from observation on BJD 2,453,101.6647 was 
rejected as it was from a twilight observation, which had both low 
precision (<trv = 78.12 m s — 1 ) and low accuracy (having a residual 
against the best Keplerian fit of over 100 m s~ 1 ). 
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Fig. 1. — Radial velocity and Keplerian model plots for the HD 37605 system. In all panels, HET ob servations are labeled with black 
filled circles, Keck observations are labeled with red crosses, and the velocities from the 2.1 m telescope (Cochran ct al. 2004) are labeled 
with blue triangles. Best Keplerian fits are plotted in black solid lines. Top left: The best-fit 2-planet Keplerian model (solid line) and the 
observed radial velocities from 3 telescopes. The HET and Keck velocities have b een a djusted to take into account the velocity offsets (i.e., 
subtracting Ajjet an d ^Kcck from the velocities, respectively; see Table [2] and § 13.21 . Bottom left: Residual velocities after subtracting 
the best-fit 2-planet Keplerian model. The lengend gives the typical size of the error bars using the ± median RV error for each telescope 
(for 2.1 m telescope only the lower half is shown). Top right: RV signal induced by HD 376056 alone, phased up to demonstrate our 
coverage. Bottom right: RV signal induced by HD 37605c alone. The two vertical dashed lines denote the date of our first observation, 
and the date when HD 37605c closes one orbit, respectively. (A color version of this figure is available in the online journal.) 

ters are calculated through bootstrapping (with 1,000 
bootstrap replicates) using the BOOTTRAN package, 
which is described in detail in the Appendi^f^l. 

The best-fit Keplerian parameters are listed in Table [2j 
The joint Keplerian fit for HD 376056 and HD 37605c has 
13 free parameters: the orbital period P, time of perias- 
tron passage T p , velocity semi-amplitude K, eccentricity 
e, and the argument of periastron referenced to the line 
of nodes tu for each planet; and for the system, the ve- 
locity offset between the center of the mass and barycen- 
ter of solar system 7 and two velocity offsets between 
the three telescopes (A Kec k and A H et, with respect to 
the velocities from the 2.1 m telescope as published in 
iCochran et"alll2004[ ). We did not include any stellar jit- 
ter or radial velocity trend in the fit (i.e., fixed to zero). 
The radial velocity signals and the best Keplerian fits 
for the system, HD 376056 only, and HD 37605c only are 
plotted in the three panels of Fig. [TJ respectively. 

Adopting a stellar mass of M* = 1.000 ± 0.017 M© (as 
in Table [I]), we estimated the minimum mass (Msini) 
for HD 376056 to be 2.802 ± 0.011 M Jup and 3.366 ± 
0.072 M Jup for HD 37605c. While HD 376056 is on a 



- 3.4 




19 The BOOTTRAN package is made publicly available online 
at |http: / / ex oplancts.or g/ code/ 1 and the Astrophysics Source Code 
Library! 



7.0 7.2 7.4 7.6 7.8 8.0 8.2 8.4 
Period (year) 

Fig. 2. — xi niap for the best Keplerian fits with fixed values 
of period P and minimum planet mass Msini for HD 37605c. 
This is showing that both P and Ms'mi are well-constrained for 
this planet. The levels of the contours mark the la (68.27%), 2a 
(95.45%) and 3a (99.73%) confidence intervals for the 2-D x 2 dis- 
tribution. 

close-in orbit at a = 0.2831 ± 0.0016 AU that is highly 
eccentric (e = 0.6767 ± 0.0019), HD 37605c is found to 
be on a nearly circular orbit (e = 0.013 ± 0.015) out at 
a = 3.814 ± 0.058 AU, which qualifies it as one of the 
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TABLE 2 
KEPLERIAN FIT PARAMETERS 



Parameter 


HD 376056 


HD 37605c 


P (days) 


55.01307 ± 0.00064 


2720 ± 57 


T p (BJD) a 


2453378.241 ± 0.020 


2454838 ± 581 


T c (BJD) b 


2455901.361 ± 0.069 




K (m s- 1 ) 


202.99 ± 0.72 


48.90 ± 0.86 


6 




n m ^ -i- n m ^ 

u.u±o ZC U.U10 


ui (deg) 


220.86 ± 0.28 


221 ± 78 


M sin i (Mj up ) 


2.802 ± 0.011 


3.366 ± 0.072 


a (AU) 


0.2831 ± 0.0016 


3.814 ± 0.058 


7 (m s~ x ) 


-50.7 ± 4.( 




A K eck (m s _1 ) c 


55.1 ± 4.7 




A H et (m s" 1 ) 


36.7 ± 4.7 






2.28 {d.o.f. = : 


L24) 


rms (m s — 1 ) 


7.61 




Jitter (m s" 1 ) 11 


3.6 





a Time of Pcriastron passage. 

Time of conjunction (mid-transit, if the system transits). 
c Offset with respect to the velocities from the 2.1 m telescope. 

If a jitter of 3.6 m s~ 1 is added in quadrature to all RV errors, xt 
becomes 1.0. 



"Jupiter analogs". 

In order to see whether the period and mass of 
the outer planet, HD 37605c, are well constrained, we 
mapped out the \t values for the best Keplerian fit in 
the P c -M c sin i space (subscript 'c' denoting parameters 
for the outer planet, HD 37605c). Each xl value on the 
P c -M c sin i grid was obtained by searching for the best-fit 
model while fixing the period P c for the outer planet and 
requiring constraints on K c and e c to maintain M sin i 
fixed. As shown in Fig. [2j our data are sufficient to have 
both P c and M c sin i well-constrained. This is also con- 
sistent with the tight sampling distributions for P c and 
M c sini found in our bootstrapping results. 

The rms values against the best Keplerian fit are 7.86 
m s" 1 for HET, 2.08 m s" 1 for Keck, and 12.85 m s" 1 for 
the 2.1 m telescope. In the case of HET and Keck, their 
rms values are slightly larger than their typical reported 
RV errors (~5m s _1 and ~ 1 m s _1 , respectively). This 
might be due to stellar jitter or underestimated system- 
atic errors in the velocities. We note that the x v is re ~ 
duced to 1.0 if we introduce a stellar jitter of 3.6 m s^ 1 
(added in quadrature to all the RV errors). 

3.3. Comparison with MCMC Results 

We compared our best Keplerian fit from RVLIN 
and uncertainties derived from BOOTTRAN (abbrevi- 
ated as RVLIN+BOOTTRAN he reafte r ) with , that from 
a Bay esian framework following Ford (2005) and iFordl 
(2006) (referred to as the MCMC analysis hereafter). Ta- 
ble |3] lists the major orbital parameters from both meth- 
ods for a direct comparison. Fig. [3] illustrates this com- 
parison, but with the MCMC results presented in terms 
of 2-D confidence contours for P, e, K, Msini, and cj of 
both planets, as well as for T c of HD 376056. 

For the Bayesian analysis, we assumed priors that 
are uniform in log of orbital period, eccentricity, argu- 
ment of pericenter, mean anomaly at epoch, and the 
velocity zero-point. For the velocity amplitude (K) 
and jitter (<7j), we adopted a prior of the form p(x) = 
(x + x a ) [log(l + x/xq)]^ 1 , with K Q — aj :Q = lm s _1 , 



i.e. high values are penalized. For a detailed discussion 
of priors, strategies to deal with correlated parameters, 
the choice of the proposal transition probability distri- 
bution function, and other details of the algorithm, w e 
refer the reader to the o riginal papers: lFoTd1|2005l[2006T) : 
iFord fe Gregory! ([2007). The likelihood for radial veloc- 
ity terms assumes that each radial velocity observation 
(iij) is independent and normally distributed about the 
true radial velocity with a variance of erf + en , where o~i is 
the published measurement uncertainty, aj is a jitter pa- 
rameter that accounts for additional scatter due to stel- 
lar variability, instrumental errors and/or inaccuracies in 
the model (i.e., neglecting planet-planet interactions or 
additional, low amplitude planet signals). 

We used an MCMC method based upon Keplerian or- 
bits to calculate a sample from the posterior distribu- 
tion (|Fordl[2006l ). We calculated 5 Markov chains, each 
with ~ 2 x 10 8 states. We discarded the first half of 
the chains and calculate Gelman-Rubin test statistics 
for each model parameter and several ancillary variables. 
We found no indications of non-convergence amongst the 
individual chains. We randomly drew 3 x 10 4 solutions 
from the second half of the Markov chains, creating a 
sample set of the converged overall posterior distribu- 
tion of solutions. We then interrogated this sample on 
a parameter-by-parameter basis to find the median and 
68.27% (lcr ) values reported in Table [3l We refer to this 
solution set below as the "best-fit" MCMC solutions. 

We note that the periods of the two planets found in 
this system are very widely separated (P c /Pb ~ 50), 
so we do not expect planet-planet interactions to be 
strong, hence we have chosen to forgo a numerically 
intensive N-bo dy DEMCMC fitting proce dure (see e.g. 
I Johnson et all 1201 lab iPavne fe Fordl |2011[ ) as the non- 
Keplerian perturbations should be tiny (detail on the 
magnitude of the perturbations is provided in ^3.4[) . 
However, to ensure that the Keplerian fits generated 
are stable, we took the results of the Keplerian MCMC 
fits and i njected those sys tems into the Mercury n-body 
package (|Chamberslll999| ) and integrated them forward 
for ~ 10 8 years. This allows us to verify that all of the 
selected best-fit systems from the Keplerian MCMC anal- 
ysis are indeed long-term stable. Further details on the 
dynamical analysis of the system can be found in i)3.4l 

We assumed that all systems are coplanar and edge-on 
for the sake of this analysis, hence all of the masses used 
in our n-body analyses are minimum masses. 

As shown in Table [3] and Fig. [3j the parameter es- 
timates from RVLIN+BOOTTRAN and MCMC meth- 
ods agree with each other very well (all within 1-cr er- 
ror bar). In some cases, the MCMC analysis reports 
error bars slightly larger than bootstrapping method 
(~ 20% at most). We note that the relatively large 
MCMC confidence intervals are not significantly reduced 
if one conducts an analysis at a fixed jitter level (e.g. 
a j = 3.5m s -1 ) unless one goes to an extremely low jitter 
value (e.g. - 1.5m s _1 ). That is, the larger MCMC er- 
ror bars do not simply result from treating the jitter as a 
free parameter. For the uncertainties on minimum planet 
mass Msini and semi-major axes a, the MCMC analy- 
sis does not incorporate the errors on the stellar mass 
estimate. Note here, as previously mentioned in § 13.11 
that the "best-fit" parameters reported by the MCMC 
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Fig. 3. — Comparison between the Bayesian (MCMC) analysis and RVLIN+BOOTTRAN results. Top four and bottom left: Contours 
of the posterior distributions of selected orbital parameters (P, e, K, M sini, and ui) based on the MCMC analysis (dashed dotted line). 
The x-axes are orbital parameters of the inner planet, b, and the y-axes are those of the outer planet, c. The inner contours mark the 68.27% 
(Ter') 2-D confidence regions and the outer ones are 95.45% ('2<r') ones. Also plotted are the best Keplerian fit from RVLIN (blue squares) 
and ±lcr error bars estimated via bootstrapping (blue bars). Bottom right: Marginalized posterior distribution of time of conjunction 
(mid-transit) T c of HD 376056 in dashed dotted line. The solid grey vertical line is the median of the distribution, and the dashed grey 
vertical lines mark la co nfide nce interval. The blue square and its error bars are for the best estimate of T c from RVLIN+BOOTTRAN 
and its ±1<t errors. See ^ I3.3I for details. (A color version of this figure is available in the online journal.) 



analysis here listed in Table [3] are not a consistent set, 
as the best estimates were evaluated on a parameter-by- 
parameter basis, taking the median from marginalized 
posterior distribution of each. Assuming no jitter, The 
best Keplerian fit from RVLIN has a reduced chi-square 
value xt — 2.28, while the MCMC parameters listed in 
Table [3] give a higher xl value of 2.91. 



3.4. Dynamical Analysis 



We used the best-fit Keplerian MCMC parameters as 
the basis for a set of long-term numerical (n-body) in- 
tegrations of the HD 37605 system u sing the Mercury 
integration package (jChambersj I1999D . We used these 
integrations to verify that the best-fit systems: (i) arc 
long-term stable; (ii) do not exhibit significant variations 
in their orbital elements on the timescale of the observa- 
tions (justifying the assumption that the planet-planet 
interactions are negligible) ; (iii) do not exhibit any other 
unusual features. We emphasize again that the planets in 
this system are well separated and we do not expect any 



TABLE 3 

COMPARISON WITH MCMC RESULTS 



Parameter 


HD 37605b 




HD 37605c 






RVLIN+BOOTTRAN 


MCMC 1 


RVLIN+BOOTTRAN 


MCMC a 



P (days) 


55.01307 ± 0.00064 


55.01250 +0.00073 -0.00075 


2720 ± 57 


2707 +57 -42 


T p (BJD) 


2453378.243 ± 0.020 


2453378.243 +0.025 -0.024 


2454838 ± 581 


2454838 +354 -435 


T c (BJD) 


2455901.361 ± 0.069 


2455901.314 +0.077 -0.081 






K (m s" 1 ) 


202.99 ± 0.72 


203.91 +0.92 -0.88 


48.90 ± 0.86 


48.93 +0.82 -0.82 


e 


0.6767 ± 0.0019 


0.6748 +0.0022 -0.0023 


0.013 ± 0.015 


0.025 +0.022 -0.017 


lj (deg) 


220.86 ± 0.28 


220.75 +0.33 -0.32 


221 ± 78 


223 +50 -52 


M (deg) b 


62.31 ± 0.15 


62.27 +0.18 -0.18 


117 ± 78 


118 +56 -51 


Msini (M Jup ) 


2.802 ± 0.011 


2.814 +0.012 -0.012 


3.366 ± 0.072 


3.348 +0.065 -0.062 


a (AU) 


0.2831 ± 0.0016 


0.2833364 +0.0000027 -0.0000027 


3.814 ± 0.058 


3.809 +0.053 -0.040 


Jitter (m s _1 ) c 


3.6 


2.70 +0.53 -0.46 







Median values of the marginalized posterior distributions and the 68.27% ('ler') confidence intervals. 
Mean anomaly of the first observation (BJD 2,453,002.671503). 

Like RVLIN, BOOTTRAN assumes no jitter or fixes jitter to a certain value, while MCMC treats it as a free parameter. See § 13.31 
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instability to occur: for the masses and eccentricities in 
question, a planet at at, ~ 0.28 AU will have c ompanion 
orbits which are Hill stable for a > 0.83 AU (|Gladmanl 



1993s), so while Hill stability does not preclude outward 
scatter of the outer planet, the fact that a c ~ 3.8 » 0.83 
AU suggests that the system will be far from any such 
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instability. 

We integrated the systems for > 10 8 years (~ 10 7 x 
the orbital period of the outer planet and > 10 2 x the 
secular period of the system), and plot in Fig. @]the evo- 
lution of the orbital elements a, e, & u>. On the left-hand 
side of the plot we provide short-term detail, illustrat- 
ing that over the ~ 10 year time period of our observa- 
tions, the change in orbital elements will be very small. 
On the right-hand side we provide a much longer-term 
view, plotting 10 7 out of > 10 8 years of system evolution, 
demonstrating that (i) the secular variation in some of 
the elements (particularly the eccentricity of the outer 
planet; see e c in red) over a time span of ~ 4 x 10 5 years 
can be significant: in this case we see 0.03 < e c < 0.11, 
but (ii) the system appears completely stable, as one 
would expect for planets with a period ratio P c /Pb ~ 50. 
Finally, at the bottom of the figure we display the range 
of parameter space covered by the et coscji, sincjj pa- 
rameters (i = b in blue for inner planet and i = c in red 
for outer planet), demonstrating that the orbital align- 
ments circulate, i.e. they do not show any signs of res- 
onant confinement, which confirms our expectation of 
minimal planet-planet interaction as mentioned before. 

As noted above, our analysis assumed coplanar plan- 
ets. As such the planetary masses used in these dynam- 
ical simulations are minimum masses. We note that for 
inclined systems, the larger planetary masses will cause 
increased planet-planet perturbations. To demonstrate 
this is still likely to be unimportant, we performed a 10 8 
year simulation of a system in which 1/ sini = 10, push- 
ing the planetary masses to ~ 30 Mj up . Even in such a 
pathological system the eccentricity oscillations are only 
increased by a factor of ~ 2 and the system remains com- 
pletely stable for the duration of the simulation. 

We also performed a separate Transit Timing Variation 
(TTV) analysis, using the best-fit MCMC systems as the 
basis for a set of highly detailed short-term integrations. 
From these we extracted the times of transit and found 
a TTV signal - 100 s, or - 0.001 day, which is much 
smaller than the error bar on T c (~ 0.07 day). Therefore 
we did not take into account the effect of TTV when 
performing our transit analysis in the next section. 

3.5. Activity Cycles and Jupiter Analogs 

The coincidence of the Solar activity cycle period of 
11 years and Jupiter's orbital period near 12 years illus- 
trates how activity cycles could, if they induced apparent 
line shifts in disk-integrated stellar spectra, confound at- 
tempts to detec t Jupi ter analogs around Sun-like stars. 
Indeed. iDravinsI (|1985l ) predicted apparent radial velocity 
variations of up to 30 m s" 1 i n sola r lines due to the So- 
lar cycle, and iDeming et al.l ()1987f ) reported a tentative 
detection of such a signal in NIR CO lines of 30 m s _1 in 
just 2 years, and noted that such an effect would severely 
hamper searches for J upiter analogs. That c oncern was 
further amplified when lCampbell et al.l (|1991l) reported a 
positive correlation between radial velocity and chromo- 
spheric activity in the active star k 1 Cet, with variations 
o f order 50-100 m s -1 . 

iWright et all (|2008l ) found that the star HD 154345 
has an apparent Jupiter analog (HD 154345 &), but that 
this star also shows activity variations in phase with 
the radial velocity variations. They noted that many 



Sun-like stars, including the precise radial velocity stan- 
dard star HD 185144 (a Dra) show similar activity vari- 
ations and that rarely, if ever, are these signals well- 
correlated with signals similar in strength to that seen in 
HD 154345 (~ 15 m s _1 ), and concluded that the simi- 
larity was therefore likely just an inevitable coincidence. 
Put s uccinctly, actiyity cy cles in Sun-like stars are com- 
mon (jBaliunas et al.lfl995l) , but few Jupiter analogs have 
been discovered, meaning that the early concern that ac- 
tivity cycles would mimic giant planets is not a severe 
problem. 

Nonetheless, there is growing evidence that activity cy- 
cles can, in some stars, induce radial velocity variations 
dSantos et al.l 120101: lLovis et all 120111: IDumusaue et all 
120121) . and the example of HD 1 54345 still warrants care 
and concern. Most significantly, IDumusaue et al.l ()2011[ ) 
found a positive correlation between chromospheric ac- 
tivity and precise radial velocity in the average measure- 
ments of a sample of HARPS stars, and provided a for- 
mula for predicting the correlation strength as a function 
of the metallicity and effective temperature of the star. 
Their formulae predict a value of 2 m s _1 for the most 
suspicious case in the literature, HD 154345 (compared to 
an actual semiamplitude of ^ 15 m s -1 ), but are rather 
uncertain. It is possible that in a few, rare cases, the 
formula might significantly underestimate the amplitude 
of the effect. 

The top panel of Fig.[3]plots the T12 APT observations 
from all five observing seasons (data provided in Table 
IU see details on APT photometry in § 14. ip . The dashed 
line marks the mean relative magnitude (A(6 + y)/2) of 
the first season. The seasonal mean brightness of the 
star increases gradually from year to year by a total of 
~ 0.002 mag, which may be due to a weak long-term 
magnetic cycle. However, no evidence is found in support 
of such a cycle i n the Mount Wilson chrom ospheric Ca II 
H & K indices (jlsaacson fc Fischer! I2010H although the 
S values vary by approxim ately 0.1 oy e r the span of a 
few years. The formulae of lLovis et all (|2Qlll ) predict a 
corresponding RV variation of less than 2 m s _1 due to 
activity, far too small to confound our planet detection 
with K = 49 m s -1 . 

Since we do not have activity measurements for this 
target over the span of the outer planet's orbit in HD 
37605, we cannot definitively rule out activity cycles as 
the origin of the effect, but the strength of the outer 
planetary signal and the lack of such signals in other 
stars known to cycle strongly dispels concerns that the 
longer signal is not planetary in origin. 

4. THE DISPOSITIVE NULL DETECTION OF TRANSITS 
OF HD 376056 

We have performed a transit search for the inner planet 
of the system, HD 376056. This planet has a tran- 
sit probability of 1.595% and a predicted transit dura- 
tion of 0.352 day, as derived from the stellar parame- 
ters listed in Table [1] and the orbital parameters given 
in Table [5J From the minimum planet mass (M sin i — 
2.802 ± 0.011 Mju D ; see Table H} and the models of 
IBodenheimer et all (|2003l ). we estimate its radius to be 
R p = 1.1 i?j up - Combined with the stellar radius of HD 
37605 listed in Table Q] R* = 0.901 ± 0.015 R Q , we esti- 
mate the transit depth to be 1.877% (for an edge-on tran- 
sit, % = 90°). We used both ground-based (APT; $0} 
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Fig. 5. — Photometric observations of HD 37605 acquired over 
five years with the T12 0.8m APT. The top panel shows the entire 
five-year data set; the dotted line represents the mean brightness of 
the first observing season. A long-term brightening trend is evident 
with a total range in the seasonal means of 0.002 mag. The middle 
panel shows the photometric data normalized so that each season 
has the same mean as the first and then phased to the orbital 
period of HD 37605b (55.01307 day). The solid line is the predicted 
transit light curve, with Phase 0.0 being the predicted time of mid- 
transit, T c . A least-squares sine fit of the phased data produces the 
very small semi-amplitude of 0.00031 ± 0.00011 mag, providing 
strong evidence that the observed radial-velocity variations are not 
produced by rotational modulation of surface activity on the star. 
The bottom panel plots the observations near T c at an expanded 
scale on the abscissa. The horizontal bar below the transit window 
represents the ±lcr uncertainty in T c . Unfortunately, none of the 
APT observations fall within the predicted transit window, so we 
a re u nable to rule out transits with the APT observations. See 
§ 14. II for more. 

and space-based (MOST; ^4.2[) facilities in our search. 

4.1. APT Observations and Analysis 

The T12 0.8-m Automatic Photoelectric Telescope 
(APT), located at Fairborn Observatory in south- 
ern Arizona, acquired 696 photometric observations of 
HP 37605 b etween 2008 January 16 and 2012 April 7. 
IHenrvl (| 1999T) provides detailed descriptions of observing 
and data reduction procedures with the APTs at Fair- 
born. The measurements reported here are differential 
magnitudes in A(6 + y)/2, the mean of the differential 
magnitudes acquired simultaneously in the Stromgren b 
and y bands with two separate EMI 9124QB bi-alkali 
photomultiplier tubes. The differential magnitudes are 
computed from the mean of three comparison stars: 
HD 39374 (V = 6.90, B-V = 0.996, KO III), HD 38145 
(V = 7.89, B-V = 0.326, FO V), and HD 38779 (V = 
7.08, B-V = 0.413, F4 IV). This improves the precision 
of each individual measurement and helps to compensate 
for any real microvariability in the comp stars. Inter- 
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Fig. 6. — Brightness variability in HD 37605 possibly induced 
by stellar rotation at P = 57.67 ± 0.30 days. Top panel is the 
periodogram of the complete, normalized data set. Bottom panel 
shows the normalized photometry folded with this possible rotation 
period. The peak-to-peak amplitude is 0.00120 ± 0.00021 mag. 
See ^ 14.11 for more. 



TABLE 4 

PHOTOMETRIC OBSERVATIONS OF 
HD 37605 FROM THE T12 0.8m APT 



Heliocentric Julian Date 
(HJD - 2,400,000) 


A(b + y)/2 
(mag) 


54,481.7133 


1.4454 


54,482.6693 


1.4474 


54,482.7561 


1.4442 


54,483.6638 


1.4452 


54,495.7764 


1.4469 


54,498.7472 


1.4470 



Note. — This table is presented in its 
entirety in the electronic edition of the As- 
trophysical Journal. A portion is shown here 
for guidance regarding its form and content. 

comparison of the differential magnitudes of these three 
comp stars demonstrates that all three are constant to 
0.002 mag or better from night to night, consistent with 
typical sing le- measurem ent precision of the APT (0.0015- 
0.002 mag: lHenr^[l999l) . 

Fig. [5] illustrates the APT photometric data and our 
transit search. As mentioned in § 13.51 the top panel 
shows all of our APT photometry covering five observ- 
ing seasons, which exhibits a small increasing trend in 
the stellar brightness. To search for the transit signal 
of HD 376056, the photometric data were normalized so 
that all five seasons had the same mean (referred to as 
the "normalized photometry" hereafter). The data were 
then phased at the orbital period of HD 376056, 55.01307 
days, and the predicted time of mid-transit, T c , defined 
as Phase 0. The normalized and phased data are plot- 
ted in the middle panel of Fig. [5] The solid line is the 
predicted transit light curve, with the predicted transit 
duration (0.352 day or 0.0064 phase unit) and transit 
depth (1.877% or ~ 0.020 mag) as estimated above. The 
scatter of the phased data from their mean is 0.00197 
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MOST photometry 



TABLE 5 

PHOTOMETRIC OBSERVATIONS OF 
HD 37605 ON MOST 



Heliocentric Julian Date 
(HJD - 2,451,545) 



Relative Magnitude 
(mag) 



4355.5105 


-0.0032 


4355.5112 


-0.0047 


4355.5119 


-0.0018 


4355.5126 


-0.0026 


4355.5133 


-0.0018 


4355.5140 


-0.0039 



Note. — This table is presented in its 
entirety in the electronic edition of the As- 
trophysical Journal. A portion is shown here 
for guidance regarding its form and content. 

mag, consistent with APT's single-measurement preci- 
sion, and thus demonstrates that the combination of our 
photometric precision and the stability of HD 37605 is 
easily sufficient to detect the transits of HD 37605b in our 
phased data set covering five years. A least-squares sine 
fit of the phased data gives a very small semi-amplitude 
of 0.00031 ± 0.00011 mag (consistent with zero) and 
so provides strong evidence that the observed radial- 
velocity variations are not produced by rotational mod- 
ulation of surface activity on the star. 

The bottom panel of Fig. [5] plots the phased data 
around the predicted time of mid-transit, T c , at an ex- 
panded scale on the abscissa. The horizontal bar below 
the transit window represents the ±ler uncertainty on 
T c (0.138 day or 0.0025 phase unit for T c 's near BJD 
2,455,901.361; see § [372]) . The light curve appears to 
be highly clustered, or binned, due to the near integral 
orbital period (P ~ 55.01 days) and consequent incom- 
plete sampling from a single observing site. Unfortu- 
nately, none of the data clusters chance to fall within the 
predicted transit window, so we are unable to rule out 
transits of HD 37605b with the APT observations. 

Periodogram analysis of the five individual observing 
seasons revealed no significant periodicity between 1 and 
100 days. This suggests that the star is inactive and the 
observed K ~ 200 m s" 1 RV signal (for HD 376056) is 
unlikely to be the result of stellar activity. 

Analysis of the complete, normalized data set, how- 
ever, suggests a week periodicity of 57.67 ± 0.30 days 
with a peak-to-peak amplitude of just 0.0012 ± 0.0002 
mag (see Fig. [6]). We tentatively identify this as the stel- 
lar rotation period. This period is consistent with the 
projected rotational velocity of usini < 1 km s _1 derived 
from our stellar analysis desc ribed in ^2.3| It is a l so con - 
sistent with the analysis of llsaacson fc Fischer! (|2010H , 
who derived a Mount Wilson chromospheric Ca II H & 
K index of S — 0.165, corresponding to log -Rhk = — 5.03. 
Together, these results imply a rotation period > 46 
days and an age of ~ 7 Gy r (see Table [J). Similarly, 
llbukivama fc Arimotol (2002) find an age of > 10 Gyr 
using isochrones along with the Hipparcos parallax and 
space motion, supporting HD 37605's low activity and 
long rotation period. 

4.2. MOST Observations and Analysis 
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Fig. 7.— Photometric observations of HD 37605 by the MOST 
satellite, which rule out the edge-on transit of HD 376056 at a > 
IOct level. The solid line is the predicted transit light curve, and the 
dashed vertical lines are the la transit window boundaries defined 
by adding cry c (0.069 day) on both sides of the predicted transit 
window (0.352-day wide). See £j |4,2l for more details. 

As noted earlier, the near-integer period of HD 376056 
makes it difficult to observe from a single longitude. The 
brightness of the target and the relatively long predicted 
transit duration creates additional challenges for ground- 
based observations. We thus observed HD 37605 during 
2011 December 5-6 (around the predicted T c at BJD 
2, 455, 901.361 as listed in Table© with the MOST (Mi- 
crovaria bility and Oscillations of Stars) satellite la unched 
in 2003 (jWalker et al.|[2Tjol [MaFthews et alj|200l in the 
Direct Imaging mode. This observing technique is simi- 
lar to ground-based CCD photometry, allowing to apply 
traditional ap erture and PSF pr ocedures for data extrac- 
tion (see e.g. iRowe et al.l [20061 . for details). Outlying 
data points caused by, e.g., cosmic rays were removed. 

MOST is orbiting with a period of ~ 101 minutes 
(14.19 cycles per day, cd _1 ), which leads to a periodic ar- 
tifact induced by the scattered light from the earthshine. 
This signal and its harmonics are further modulated with 
a frequency of 1 cd _1 originating from the changing 
albedo of the earth. To correct for this phenomenon, 
we constructed a cubic fit between the mean background 
and the stellar flux, which was then subtracted from the 
data. The reduced and calibrated MOST photometric 
data are listed in Table [5] 

The MOST photometry is shown in Figure [7] for the 
transit window observations. The vertical dashed lines 
indicate the beginning and end of the 1 a transit window 
defined by adding <jt c (0.069 day) on both sides of the 
predicted transit duration of 0.352 days. The solid line 
shows the predicted transit model for the previously de- 
scribed planetary parameters. The mean of the relative 
photometry is 0.00% (or 0.00 mag), with a rms scatter 
of 0.17%, and within the predicted transit window there 
are 58 MOST observations. Therefore, the standard er- 
ror on the mean relative photometry is 0.17%/y / 58 = 
0.022%. This means that, for the predicted transit win- 
dow and a predicted depth of 1.877%, we can conclude 
a null detection of HD 376056's transit with extremely 
high confidence (149tr). 

Note that the above significance is for an edge-on tran- 
sit with an impact parameter of b — 0.0. A planetary 
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trajectory across the stellar disk with a higher impact 
parameter will produce a shorter transit duration. How- 
ever, the gap between each cluster of MOST measure- 
ments is 0.06 days which is 17% of the edge-on tran- 
sit duration. In order for the duration to be fit within 
the data gaps, the impact parameter would need to be 
b > 0.996. To estimate a more conservative lower limit 
for 6, we now assume the most unfortunate case where 
the transit center falls exactly in the middle of one of the 
measurement gaps, and also consider the effect of limb 
da rkening by using the n on-linear limb darkening model 
bv lMandel fc Agoll (|2002D with their fitted coefficients for 
HD 209458. Even under this scenario, we can still con- 
clude the null detection for any transit with b < 0.951 at 
> 5(7 (taking into account that there are at least ~ 20 ob- 
servations will fall within the transit window in this case, 
though only catching the shallower parts of the transit 
light curve). 

All of the above is based on the assumption that the 
planet has the predicted radius of 1.1 i?j up - If in reality 
the planet is so small that even a b = transit would fall 
below our detection threshold, it would mean that the 
planet has a radius of < 0.36 i?j up (a density of > 74.50 
g/cm 3 ), which seems unlikely. It is also very unlikely 
that our MOST photometry has missed the transit win- 
dow completely due to an ill-predicted T c . In the sam- 
pling distribution of T c from BOOTTRAN (with 1000 
replicates; see S 13.21 and Appendix), there is no T c that 
would put the transit window completely off the MOST 
coverage. In the marginalized posterior distribution of 
T c calculated via MCMC (see §013] and Fig.©, there is 
only 1 such T c out of 3 x 10 4 (0.003%). 

5. SUMMARY AND CONCLUSION 

In this paper, we report the discovery of HD 37605c 
and the dispositive null detection of non-grazing tran- 
sits of HD 376056, the first planet discovered by HET. 
HD 37605c is the outer planet of the system with a pe- 
riod of ~ 7.5 years on a nearly circular orbit (e = 0.013) 
at a = 3.814 AU. It is a "Jupiter analog" with Msini = 
3.366 -Mj up , which adds one more sample to the currently 
still small inventory of such planets (only 10 including 
HD 37605c; see §0}. The discovery and characterization 
of "Jupiter analogs" will help understanding the forma- 
tion of gas giants as well as the frequency of true solar 
system analogs. This discovery is a testimony to the 
power of continued observation of planet-bearing stars. 

Using our RV data with nearly 8-year long baseline, we 
refined the orbital parameters and transit ephemerides 
of HD 376056. The uncertainty on the predicted mid- 
transit time was constrained down to 0.069 day (at and 
near T c = 2,455,901.361 in BJD), which is small com- 
pared to the transit duration (0.352 day). In fact, just 
the inclusion of the two most recent points in our RV 
data have reduced the uncertainty on T c by over 10%. 
We have performed transit search with APT and the 
MOST satellite. Because of the near-integer period of 
HD 376056 and the longitude of Fairborn Observatory, 
the APT photometry was unable to cover the transit win- 
dow. However, its excellent photometric precision over 
five observing seasons enabled us to rule out the possi- 
bility of the RV signal being induced by stellar activity. 
The MOST photometric data, on the other hand, were 
able to rule out an edge-on transit with a predicted depth 



TABLE 6 

Updated Msini and Errors for HD 1147626 
and HD 1684436, c 



Planet 


Msini ± std. error (Mj up ) 


HD 1147626 a 


11.086 ±0.067 


HD 1147626 b 


11.069 ±0.063 


HD 1684436 


7.696 ± 0.015 


HD 168443c 


17.378 ± 0.044 



a For best orbital fit with RV trend (dv/dt). 
b For best orbital fit without RV trend 
{dv/dt). 



of 1.877% at a ^ 10er level, with a 5<r lower limit on the 
impact parameter of 6 < 0.951. This transit exclusion is 
a further demonstration of the TERMS strategy, where 
follow-up RV observations help to reduce the uncertainty 
on transit timing and enable transit searches. 

Our best-fit orbital parameters and errors from 
RVLIN+BOOTTRAN were found to be consistent with 
those derived from a Bayesian analysis using MCMC. 
Based on the best-fit MCMC systems, we performed dy- 
namic and TTV analysis on the HD 37605 system. Dy- 
namic analysis shows no sign of orbital resonance and 
very minimal planet-planet interaction. We derived a 
TTV of ^100 s, which is much smaller than ax c - 

We have also performed a stellar analysis on HD 37605, 
which shows that it is a metal rich star ([Fe/H] = 0.336± 
0.030) with a stellar mass of M* = 1.000 ± 0.017 M 
with a radius of i?* = 0.901 ± 0.015. The small variation 
seen in our photometric data (amplitude < 0.003 mag 
over the course of four years) suggests that HD 37605 is 
consistent as being an old, inactive star that is probably 
slowly rotating. We tentatively propose that the rotation 
period of the star is 57.67 ± 0.30 days, based on a weak 
periodic signal seen in our APT photometry. 

6. NOTE ON PREVIOUSLY PUBLISHED ORBITAL FITS 

In early 2012, we repaired a minor bug in the BOOT- 
TRAN package, mostly involving the calculation and er- 
ror bar estimation of Msini. As a result, the Msini 
values and their errors for two previously published sys- 
tems (three planets) need to b e updated. They are: 
HD 1147 626 (Kane et al.|[20lTal) . HD 1684436, and HD 
168443c (|Pilvavskv et aflboilbl ). Table [6] lists the up- 
dated Msini and error bars. 

One additional system, HD 63454 dKane et al.ll2011dlh 
was also analyzed using BOOTTRAN. However, the 
mass of HD 634546 is small enough compared to its host 
mass and thus was not affected by this change. 
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APPENDIX 
UNCERTAINTIES VIA BOOTSTRAPPING 

The unc ertainties listed for the orbital paramete r estimates^] and transit mid-time T c are calculated via boot- 
strapping (Frcc dma n 1981; Davison fc HinklevifTOfll using the package BOOTTRAN, which we have made publicly 
available (see § 13. 2[) . It is designed t o calculate error bars for transit ephemerides and the Keplerian orbital fit param- 
eters output by the RVLIN package (Wr ight fc Howard! 120091 ) . but can also be a stand-alone package. Thanks to the 
simple concept of bootstrapping, it is computationally very time-efficient and easy to use. 

The basic idea of bootstrap is to resample based on original data to create bootstrap samples (multiple data repli- 
cates); then for each bootstrap sample, derive orbital parameters or transit parameters through orbital fitting and 
calculation. The ensemble of parameters obtained in this way yields the approximate sampling distribution for each 
estimated parameter. The standard deviation of this sampling distribution is the standard error for the estimate. 

We caution the readers here that there are regimes in which the "approximate sampling distribution" (a frequentist's 
concept) is not an estimate of the posterior probability distribution (a Bayesian concept), and there are regimes (e.g., 
when limited sampling affects the shape of the x 2 surface) where there are qualitative differences and the bootstrap 
method dramatically underestim ates uncertainties (e.g., long-per iod planets when the observations are not yet sufficient 
to pin down the orbital period; iFordl 120051 : 1 Bender et al.l I2012D . In situations with sufficient RV data, good phase 
coverage, a sufficient time span of observations and a good orbital fit, bootstrap often gives a useful estimate of the 
parameter uncertainties. For the data considered in this paper, it was not obvious that the bootstrap uncertainty 
estimate would be accurate, as the time span of observations is only slightly longer than the orbital period of planet 
c. Nevertheless, we find good agreement between the uncertainty estimates derived from bootstrap and MCMC 
calculations. 

The radial velocity data are denoted as {t, v, a}, where each ij, Vi, ai represents radial velocity Vi observed at 
time (BJD) ti with velocity uncertainty crj. Extreme outliers should be rejected in order to preserve the validity of 
our bootstrap algorithm. We first derive our estimates for the tru e orbital parameters from the original RV data via 
orbital fitting, using the RVLIN package (|Wright fc Howard! 120011 : 

P = n(t, v, a), (1) 

where /3 is the best fitted orbital parameter^. From /3, we derive {t, Vb est (/3)} 7 the best-fit model (here t are treated 
as predictors and thus fixed). Then we can begin resampling to create bootstrap samples. 

Our resampling plan is model-based resampling, where we draw from the residuals against the best-fit model. For 
data that come from the same instrument or telescope, in which case no instrumental offset needs to be taken into 
account, we simply draw from all residuals, {v — Vbest}, with equal probability for each (vi — Vb es t,i)- This new ensemble 
of residuals, denoted as r*, is then added to the best-fit model Vbest to create one bootstrap sample, v* E3- Associated 
with r* , the uncertainties a are also re-assigned to v* - that is, if Vj — Vb es t,j is drawn as rk and added to Vk to generate 
then the uncertainty for vf, is set to be <Xj. 

For data that come from multiple instruments or multiple telescopes, we incorporate our model-based resampling plan 
to include stratified sampling. In this case, although data from each instrument or telescope are close to homoscedastic, 
the entire set of data are usually highly heteroscedastic due to stratification in instrument/telescope radial velocity 
precision. Therefore, the resampling process is done by breaking down the data into different groups, {v{, v^, ■ ■ ■}, 
according to instrument and/or telescope, and then resample within each subgroup of data with the algorithm described 
in last paragraph. The bootstrap sample is then v* = {v*i, v*2, • • •}• 

To construct the approximate sampling distribution of the orbital parameter estimates /3, we compute 

f3* = n(t, v*, a*) (2) 

for each bootstrap sample, {<, v*, er*}. The sampling distribution for each orbital parameter estimate /3.; can 

(l) (2) 

be constructed from the multiple sets of p* calculated from multiple bootstrap samples (p* , p* , . . . from 



21) 



Through out the paper and sometimes in this Appendix, we 21 A „ Hescrihed in iT31 this includes the P T K e and u, for 

c , ,i a (■ i r ,1 . ,, / t i- -lie rt_S QCSCllUCQ 111 < l O . Z I L111S 111C1L1QCS HID 1 . ± 71 , 1\ . C. cillLl LU IOI 

refer to the estimates of the parameters (as distinguished from 



each planet, as well as 7, dt>/dt (if applicable), and velocity offsets 



the "true parameters", which arc not known and can only be esti- , . ■ ', , /, , ' Vr ,■ , , *%' c , 

. „ J between instruments/telescopes (it applicable) tor the system 



mated) simply as the "parameters 



22 We simply use the raw residual instead of any form of modified 
residual, because the RV data for any single instrument or telescope 
are usually close enough to homoscedasticity. 
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v* , v* ,...)■ The standard errors for /3 are simply the standard deviations of the sampling distribution^ 2 "^!. 

The sampling distribution of the estimated transit mid-time, T c , is calculated likewise. Here T c is the transit time 
for a certain planet of interest in the system, and is usually specified to be the first transit after a designated time 
T. However, the situation is complicated by the periodic nature of T c . Our approach is to first calculate, based on 
the original RV data, T c q, the estimated mid-time of the first transit after time To (an arbitrary time within the RV 
observation time window of [min(i), max(i)]; T c o is also within this window). Then 

T C = N-P + T c0 , (3) 

where P is the best-estimated period for this planet of interest, and N is the smallest integer that is larger than 
(T — T c q)/P. Next we compute T* for each bootstrap sample {t, v* , a*}. Given that within the time window of 
radial velocity observations ([min(t), max(t)]), the phase of the planet should be known well enough, it is fair to 
assume that T c q is an unbiased estimator of the true transit mid-time. Therefore we assert that T* has to be well 
constrained and within the range of [T c q — P*/2, T c o + P*/2], where P* is the period estimated from this bootstrap 
sample. If not, then we subtract or add multiple P*'s until T* falls within the range. Then naturally 

T* = N ■ P* + T* . (4) 

The ensemble of T*'s gives the sampling distribution of T c and its standard error. Note that T* is not necessarily 
within the rage of [T c - P/2, T c + P/2]. 

Provided with the stellar mass and its uncertainty, we calculate, for each planet in the system, the standard 
errors for the semi-major axis a and the minimum mass of the planet Mp im ; n (denoted as Msmi in the main text 
as commonly seen in literature, but this is a somewhat imprecise notation). As the first step, the mass function is 
calculated for the best-fit (3 and each bootstrap sample /?*, 

HP Kr)- PK " {1 ~ 6)3/2 - (Mp-^O 3 (5) 

The sampling distribution of /(P, K, e) then gives the standard error of the mass function. The minimum mass of 
the planet Mp. m ; n is then calculated by assuming sini = 1 and solving for M p . Standard error of Mp im ; n is derived 
through simple propagation of error, as the covariance between M* and /(P, K, e) is probably negligible. 
For the semi-major axis a, 

3 _ P 2 G(M* + M p ) P 2 G(M, + M p , min ) 
a - 4^2 ~ 4^2 • W 

The standard error of P 2 is calculated from its bootstrap sampling distribution, and via simple propagation of error 
we obtain the standard error of a (neglecting covariance between P 2 , M p , m i n , and M*). 



23 The standard deviation of a sampling distribution is estimated 
in a robust way using the IDL function robustsigma, which is writ- 
ten by H. F ruedenreich ba s ed on the principles of robust estimation 
outlined in Hoaglin ct al. (1983). 
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TABLE 7 

HET and KECK RADIAL VELOCITIES FOR HD 37605 





Velocity 


Uncertainty 




BJD-2440000 


(m s _1 ) 


(m s _1 ) 


Telescope 


i onnn c *7i cno 

loUUz.o715Uo 


122.4 


8.8 


HET 


i onno coco/17 

loUUo.o85247 


126.9 


5.6 


HET 


13006.662040 


132.1 


5.2 


HET 


i onno t* i n tr n 

13008.664059 


130.4 


4.8 


HET 


i om n on a '~?r-'~7 

13010.804767 


113.4 


4.7 


HET 


13013.793987 


106.2 


4.9 


HET 


"1 'Willi i \ ■-- ( \ / • - I 

13042.727963 


-120.0 


4.7 


HET 


i on/ -1 "i tr i 

13061.667551 


126.5 


4.0 


HET 


13065.646834 


111.3 


5.8 


HET 


i on^i c/noon 

13071.643820 


106.7 


5.5 


HET 


13073.638180 


91.7 


4.6 


HET 


i onon cooTnn 

1oU8z.o2o7U9 


53.4 


5.7 


HET 


lOflOQ K O C. Q C. 7 

loUoo.oyooo * 


45. 1 


6.4 


HrL 1 


i onoo cno f 7 , 7c: 

13088.593775 


30.7 


11.3 


HET 


i onon rficrcn 

13089.595750 


0.8 


5.5 


HET 


i onnn crn r 7noo 

13092.597983 


-25.5 


6.2 


HET 


i onn ,1 cro/^CTTn 

13094.586570 


-51.0 


5.9 


HET 


i onn tr cor /i 1 o 

1oU95.58o41o 


-72.2 


6.2 


HET 


loUyo.o8f4oz 


-75.0 


9.3 


rLrjj 1 


i onno t'7/'ni o 

13098.576213 


-173.1 


6.4 


HET 


i o nr ,1 nc 1 ore 

13264.951365 


-194.0 


9.8 


HET 


1 oorc n /I ""7 /I OO 

13265.947438 


-238.2 


10.5 


HET 


13266.945977 


-264.2 


13.1 


HET 


13266.959470 


-286.8 


11.6 


HET 


1o28o.9z24U7 


118.8 


9.6 


HET 


i ooi o 01 nnm 

13318.819260 


-142.8 


6.6 


HET 


i o o o p" nni' ta 7 p rn 

13335.921770 


143.4 


6.3 


HET 


i oooo nnmno 

13338.906008 


132.1 


5.4 


HET 


TOOT 1 ? 01 A/I AC 

loo77.8194U5 


-279.6 


6.1 


HET 


1 Q97Q Ol 1 oon 

loo ro.ollooU 


-161.3 


5.2 




i o o 7n onoo H 7 

13379.802247 


-39.7 


5.2 


HET 


1 O O O "1 r" A A <~>0 A 

Ioo81.o44284 


74.4 


4.7 


HET 


TOOO/I r/irCTOO 

13384.646538 


112.5 


5.8 


HET 


1 O 70 ,1 OCCOOI 

13724.855831 


88.1 


5.3 


HET 


13731.697227 


57.8 


5.0 


HET 


i otoo /' f 7/i r 7nn 

13738.674709 


36.4 


4.9 


HET 


1 otjo oTmnr 

13743.810196 


14.1 


5.7 


HET 


Io748.o472o4 


-17.8 


5.6 


HET 


14039.850147 


-104.1 


5.5 


HET 


14U54. yo45oo 


68.8 


7.3 


HET 


14055.952778 


49.3 


6.5 


HET 


14067.762810 


22.1 


6.0 


HET 


i zioT/i n n i n o i • 

14o74.924U8u 


32.7 


8.6 


HET 


14o94.8o44o8 


-12.1 


8.5 


HET 


aaa /in ooonno 

1444U.88o298 


11.9 


8.5 


HET 


144o7.82o5oz 


-121.9 


6.3 


HET 


1448 / .oloUoz 


40.4 


5.6 


HrL 1 


14515. o89o41 


-55.6 


5.7 


HET 


i a cr crn mnom 

14550.600301 


23.1 


5.5 


HET 


14759.878895 


34.3 


4.9 


HET 


1 J77C OOm 07 

14776.839187 


2.0 


6.7 


HET 


14787. 79oo98 


-40.8 


5.1 


HET 


1 4QD7 (\A 7^99 


-94.4 


5.1 


HET 


15089.968860 


76.4 


4.1 


HET 


15104.917965 


54.0 


4.6 


HET 


15123.882222 


-44.4 


4.6 


HET 


15172.734743 


-0.4 


4.6 


HET 


15172.749345 


3.0 


4.9 


HET 


15180.715307 


-65.0 


5.3 


HET 


15190.692697 


-257.0 


4.9 


HET 


15190.704818 


-252.3 


6.3 


HET 
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TABLE 7 — Continued 



Velocity Uncertainty 
BJD-2440000 (m s" 1 ) (m s" 1 ) Telescope 



15195.678923 


-47.7 


5.1 


HET 


15200.809909 


85.1 


4.5 


HET 


1 run/ 1 /~' r' f c 1 c 

15206.665616 


75.9 


4.5 


HET 


15208.646672 


72.6 


5.7 


HET 


15266.634586 


81.1 


4.6 


HET 


15274.598361 


57.1 


4.6 


HET 


15280.605398 


24.3 


5.6 


HET 


15450.974863 


25.2 


5.2 


HET 


15469.928975 


-81.7 


4.4 


HET 


15481.884352 


110.8 


4.8 


HET 


15494.004309 


71.2 


4.9 


HET 


15500.853302 


41.1 


6.0 


HET 


15505.824457 


21.2 


4.8 


HET 


15510.809500 


-13.9 


5.5 


HET 


15511.813809 


-18.4 


4.7 


HET 


15518.798148 


-139.1 


4.9 


HET 


15518.943518 


-148.8 


4.6 


HET 


15524.791579 


-97.7 


6.0 


HET 


15526.759230 


51.2 


5.8 


HET 


15526.767959 


56.5 


6.1 


HET 


15527.921902 


89.2 


4.6 


HET 


15528.921753 


105.1 


5.5 


HET 


15530.757847 


120.6 


5.3 


HET 


15532.760658 


117.6 


4.8 


HET 


15544.730827 


81.8 


4.9 


HET 


15545.871710 


87.7 


5.2 


HET 


15547.717980 


78.6 


6.5 


HET 


15547.720596 


81.2 


7.8 


HET 


15550.703130 


75.8 


5.5 


HET 


15556.688398 


43.4 


5.0 


HET 


15557.842853 


40.1 


4.9 


HET 


15566.657803 


-18.4 


5.6 


HET 


15566.666973 


-19.0 


4.9 


HET 


15569.670948 


-50.7 


4.5 


HET 


15571.659560 


-94.0 


7.4 


HET 


15571.662188 


-99.5 


7.5 


HET 


15582.777052 


75.1 


5.4 


HET 


13982.116400 


-312.6 


1.1 


Keck 


13983.110185 


-291.6 


1.1 


Keck 


13984.104595 


-175.0 


1.1 


Keck 


14024.076794 


-39.9 


1.2 


Keck 


14129.848461 


-17.1 


1.1 


Keck 


14131.926354 


-32.2 


1.3 


Keck 


14138.913472 


-98.8 


1.2 


Keck 


14544.815162 


54.2 


1.0 


Keck 


14545.799155 


50.4 


1.0 


Keck 


14546.802812 


48.1 


1.0 


Keck 


15172.939733 


21.7 


1.1 


Keck 


"1 P" .1 O P" 1 O O C f \ T* 

15435.133695 


97.5 


1.0 


Keck 


15521.851677 


-244.2 


1.2 


Keck 


15522.947381 


-268.1 


1.1 


Keck 


loozo.yooiuo 


1 O.D 


i i 
i . ± 


1YCCK 


15528.877669 


118.4 


1.0 


Keck 


15543.032113 


115.5 


1.2 


Keck 


15545.105481 


110.3 


1.2 


Keck 


15546.046328 


110.6 


1.2 


Keck 


15556.071287 


69.6 


1.1 


Keck 


15556.954989 


66.4 


1.1 


Keck 


15558.027806 


62.3 


1.2 


Keck 


15584.769501 


128.2 


1.2 


Keck 


15606.952224 


91.6 


1.2 


Keck 


15634.789262 


-81.0 


1.0 


Keck 
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TABLE 7 



Continued 



BJD-2440000 



Velocity Uncertainty 

(m s" 1 ) (m s _1 ) Telescope 



15636. 
15636. 
15663. 
15671. 
15672. 
15673. 
15793. 
15843. 



793990 
799360 
738831 
731754 
732058 
745070 
124325 
017282 



80.7 
85.6 
87.6 
45.8 
37.1 
27.6 
-86.7 
-6.4 



1.0 
1.1 
1.0 
1.2 
1.0 
1.0 
1.2 
1.3 



Keck 
Keck 
Keck 
Keck 
Keck 
Keck 
Keck 
Keck 



